Computer simulation of quantum melting in hydrogen clusters 
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We introduce a new criterion — based on multipole dynamical correlations calculated within Reptation Quan- 
tum Monte Carlo — to discriminate between a melting vs. freezing behavior in quantum clusters. This criterion 
is applied to small clusters of para-hydrogen molecules (both pristine and doped with a CO cromophore), for 
cluster sizes around 12 molecules. This is a magic size at which para-hydrogen clusters display an icosahedral 
structure and a large stability. In spite of the similar geometric structure of CO@(pH2)i2 and (pH2)i3, the 
first system has a rigid, crystalline, behavior, while the second behaves more like a superfluid (or, possibly, a 
super solid). 
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Understanding the dynamics of quantum many-body sys- 
tems is one of the major challenges presently set to theoretical 
and computational condensed-matter physicists. The combi- 
nation of density-functional theory with molecular dynamics 
realized 20 years ago by Car and Parrinellc 1 opened the way 
to the study of the dynamics of quantum driven classical sys- 
tems (i.e. of atomic systems whose dynamics is essentially 
classical, but driven by quantum-mechanical forces). The de- 
velopment of density-functional perturbation theory has al- 
lowed for a systematic calculation of the low-lying quan- 
tum excited states of these same systems in the harmonic 
approximation 2 . Quantum Monte Carlo methods, on the other 
hand, have been very successful in describing the ground-state 
and finite-temperature properties of interacting bosons 3 and 
of lattice models of strongly interacting fermionsi, and they 
promise a similar success in the study of chemical systems 
in the near future 5 . In spite of all these progresses, the abil- 
ity to calculate in a reliable way the properties of the excited 
states of interacting quantum systems remains a largely un- 
achieved goal. Recent advances in the quantum Monte Carlo 
methodology have partially changed this scenario, at least in 
what concerns systems whose ground-state wave-function is 
positive (i.e. bosons) and whose low-lying energy spectrum is 
dominated by few excited states, which is the typical situation 
for superfluids 6,7 . Thanks to these advances, it is now possi- 
ble to calculate the low-lying excitation spectrum of clusters 
of up to a few tens 4 He atoms, possibly doped with some cro- 
mophore molecules which are experimentally used as spec- 
troscopic probes of the dynamical and superfluid properties of 
the droplet. 

Molecular hydrogen, in its nuclear-spin realization called 
para-hydrogen (pH2, / = J = 0), is the only substance occur- 
ring in nature, other than 4 He, which can possibly exhibit the 
phenomenon of superfluidity at (not too) low temperature 8 . 
In spite of the lighter mass with respect to 4 He, the inter- 
molecular potential is so much stronger that the estimated 
value of the A-transition temperature 2°K) is much lower 
than the observed triple-point temperature (13.96°K). For this 
reason, much attention has been and is being paid to those 



confined geometries (such as clusters^iiSiiiii 2 ^ or filmsi 4 *!^) 
which may hinder crystallization and make superfluidity (or 
what remains of it in these geometries) observable in pH2- 
Discriminating between a melting vs. freezing behavior in a 
finite system is a subtle issue which requires a careful con- 
sideration of the dynamical evolution of the system. In this 
paper we propose a strategy to cope with this problem, based 
on Reptation Quantum Monte Carlo (RQMC), a method that 
we developed a while ago to deal with dynamical correlations 
in systems of interacting bosons 6 . As a (still very prelimi- 
nary) application, we show how using suitable (imaginary-) 
time correlation functions it is easy (and fun!) to discrimi- 
nate between the melting vs. freezing behavior of small 
clusters — with or without a CO molecule solvated in them — 
for sizes close to n = 13, a magic number which favors crys- 
tallization. 



Reptation Quantum Monte Carlo 

The dynamics of classical stochastic processes bears close 
similarities with the time evolution of quantum systems in 
imaginary time. These similarities — together with the fact 
that every system tends towards the ground state at large imag- 
inary times — lay at the basis of many quantum Monte Carlo 
methods which sample ground state properties from the equi- 
librium properties of a suitably defined classical random walk. 
Comparatively minor attention has been paid to the dynamical 
properties of the random walks used in quantum simulations. 
In this section we show how a careful exploitation of these 
similarities can be used to estimate the dynamical properties 
of strongly interacting boson fluids. 

Most ground-state Monte Carlo techniques are based on 
the property that the imaginary-time evolution of (almost) any 
trial wave-function, <3>o> in the infinite-time limit tends to the 
ground state, SPq. As a consequence, the ground- state energy 
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can be expressed as: 
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where H = —-§^2 + V{x) is the Hamiltonian of the system 
and Zq = (§o\e~ Hr \$o)- If %o is thought as the partition 
function of an auxiliary, fictitious, system, then Eqs. and 
© express the free and internal energies of this system in the 
zero-temperature limit. By using Eq. 0, ground-state expec- 
tation values of static operators and generalized susceptibili- 
ties can be expressed as first and second logarithmic deriva- 
tives of Zq with respect to the strength of a suitably chosen 
perturbation to the Hamiltonian^. 

A link with the theory of stochastic processes can be estab- 
lished by splitting the Hamiltonian as: 



H = H + £(x), 
Mx), £{x) = 
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where H + i^y-ov^7> - 

V(x), and the double prime indicates the Laplacian. Note that 
by construction $0 is an eigenstate of TL with zero eigenvaue 
and that, if it is chosen to be node-less, it also has to be its 
non-degenerate ground state, so that the excited-state spec- 
trum of H is strictly positive. £(x) is easily recognized as the 
definition of the local energy, well known to the variational 
Monte Carlo (VMC) and diffusion Monte Carlo (DMC) prac- 
titioners. The pseudo-partition function, Zo, can be given a 
path-integral representation by Trotter- splicing the propagator 
appearing therein: 

Zo = J $o(x N )g € (x N ,x N -i)g € (x N -i,x N -2)-- 

g e (x u x )Mxo)e- e ^ £iXn) V[X}. (5) 

The short- time propagator, Q e , can be expressed in terms of 
the transition matrix of a suitably defined Markovian random 
walk: 
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is the transition probability 
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associated with the Langevin equation: 

dx = f(x)e + d£, 

= and £( r ) is a Wiener process: (df) = 

°; (( d 2 ) = 2e - p o(x) = $o(x) 2 is the (unique) equilib- 
rium distribution of Eq. Q. This fact, together with Eq. ©, 
allows to cast Eq. Q into the form: 
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where S[X] = e £(x n ) is the reduced action of the sys- 
tem, and the brackets, (-)rw, indicate a statistical average 
over an equilibrium distribution of quantum paths (or random 
walks), {X}. 

Suppose now that the system is coupled to a (imaginary-) 
time dependent perturbation: H\ = H + A(r)A The second 
derivatives of the pseudo-partition function, Zo, with respect 
to A calculated at A = give the ground-state imaginary -time 
correlations of the A operator: 
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According to the above equations, a practical algorithm to 
calculate ground- state time correlations would be to generate 
segments of a random walk according to the Langevin equa- 
tion, Eq. Q; each segment is then accepted or rejected ac- 
cording to a Metropolis test performed on the reduced action, 
S[X]; the statistical average of time correlations calculated 
along the paths thus generated will provide the desired quan- 
tum time correlations. This is the heart of the Reptation Quan- 
tum Monte Carlo method 6 . This algorithm can be generalized 
by generating trial moves for the paths according to any a- 
priori probability distribution, and then accepting or rejecting 
them according to a Metropolis test done one the actual prob- 
ability distribution for the paths, Eq. 0, 16 . 

The time correlation function of various observables which 
couple to external fields (such as, e.g., the electric dipole) are 
the Laplace transforms of spectral weight functions which, 
in some favorable cases, can be directly compared with 
experiments 7 . In the next section we show how suitably de- 
fined time-correlation functions can be used to discriminate 
the quantum melting vs. freezing behavior of small pH2 clus- 
ters. 



Magic pH2 clusters 

Infrared and microwave spectroscopies of solvated cro- 
mophores are being widely used to probe the microscopic 
structure of 4 He nano-droplets, as well as their propensity to 
display a superfluid behavior 17 . It is hoped that the use of 
these techniques for pH2 clusters of different size will help 
identify fingerprints of superfluidity, and to clarify the inter- 
play between this phenomenon and the competing tendency 
to crystallization in these systems. Infrared spectra spectra of 
carbon mono-oxide (CO) molecules solvated in small (pH2) n 
clusters (n < 14) have been recently measured and RQMC 
simulation have helped unravel them and assign the identity 
of the individual lineA 

In Fig. 1 we display the radial distribution functions of pH2 
around a CO molecule, for cluster sizes around the magic 
number n = 13i£. For n < 12 the magnitude of the max- 
imum of these functions increases with the cluster size. For 
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FIG. 1: Left panel: radial density distribution function of pFh 
molecules in CO@(pH 2 ) n (dotted line: n — 11, continuous line: 
n — 12; dash-dotted line: n — 13). Right panel: pFb chemical 
potential (/i/v = K - K-i) as a function of the cluster size in 
CO@(pH 2 )n. 



n > 12 this value stays nearly constant and the increase of 
the cluster size shows in the tail of the distribution, rather than 
in the height of the peak, indicating that the first solvation 
shell is completed at this magic number. For n = 12 the ab- 
solute value of chemical potential also displays a maximum, 
suggesting the greatest stability of the clusters at this size. 
At the same time, the rotational distortion constant displays 
a minimum, indicating the largest rigidity of this clusters. All 
these findings suggest that the propensity of the droplet 
to crystallize is maximum ad this cluster size. This fact has 
to be compared with the findings of Ref. 9 which show clear 
signs of superfluidity in a cluster of 13 pF^ molecules, whose 
structure is similar to that of CO@(pH2) n with the solvated 
CO molecule substituted with a molecule. In order to as- 
sess the melting vs. freezing behavior of finite systems, we 
introduce a dynamical criterion based on the persistence of 
geometrical signatures of the crystal structure. 

The shape of a rigid body is well described by the multipole 
moments of, say, its mass density distribution around the cen- 
ter of mass, Q l m . In order to characterize the degree of rigidity 
of a cluster which, in addition to shape fluctuations, also ex- 
periences rotational diffusion, the shape must be described in 
a way which is independent of the orientation. The magni- 
tude of the multipole, defined as q l = Y^m Qln Qm> provides 
such an invariant characterization of a rigid body. This quan- 
tity, however, is not suitable to discriminate the cases where 
a given multipole vanishes on the average, and (q l ) ^ be- 
cause of the fluctuations, from those where a non vanishing 
value of {q l ) is due to the average shape of the cluster. This 
discrimination is best achieved in terms of a suitably defined 
rotationally-invariant time correlation functions of the various 
multipole moments in the rotating-axes frame. In order to 
proceed, let us define an effective angular velocity, fi(r), as 
the average of the forward angular velocities of the individual 
molecules. The operator corresponding to the global rotation 
accomplished in a time step e reads: 
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mentum). The total rotation accomplished during a time r is: 

(11) 
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R(t) = Te- l $>- n ^ dT ' 

T/e 

= l[AR(ne), 



where Te^ indicates the time-ordered exponential. With the 
aid R(r) we define a rotating-axes multipole time-correlation 
function, q(t), as the correlation between the multipole mo- 
ment at (imaginary) time r' and the multipole at time r' + r, 
rotated by ^(V) -1 (time-translation invariance implies that 
these correlations are independent of r')\ 
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where Q l m (r) is the multipole moment of the system at (imag- 
inary) time r, the brackets, (•), indicate an average over 
the random walk, and Q indicates the multipole rotated by 
i? -1 (r). If a system is rigid, q(t) is independent of time, 
whereas it tends to zero at large values of r whenever the 
shape of the system undergoes random fluctuations. A large 
value of ci (r) at large times is thus an indicator of the persis- 
tence of the shape of the system 19 . 

In the case of CO@(pH2)i2 our simulations indicate a very 
rapid decay of the correlations for all the multipoles up to and 
including I = 5. For I = 6, instead, after a short transient in 
which a very steep drop occurs, the correlations decay very 
slowly, indicating a long persistence of a shape compatible 
with a nonvanishing I = 6 multipole and vanishing multi- 
poles with smaller values of I. In Fig. 2 we report the I = 6 
rotating-axes multipole time correlation function calculated 
for CO@fcH 9 )i9. 
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where L is the generator of the rotation group (angular mo- 



FIG. 2: / = 6 rotating-axes multipole time correlation function (Eq. 
[13} for CO@(pH 2 )i2. 

This behavior is indeed compatible with a regular icosa- 
hedral shape, whose lowest non-vanishing multipoles corre- 
spond to I = 6 and I = 10. In order to visualize the intrin- 
sic shape of the cluster, we calculate the rotating-axes num- 
ber density distribution of the hydrogen molecule, nn 2 ( r )> by 



4 




FIG. 3: pH2 rotating- axes number- density distribution for 
CO@(pH2)i2 (see text). The displayed iso-surface corresponds to 
a value of the density including 60% of the total number of H2 
molecules. The solvated CO molecule rotates nearly freely with re- 
spect to the rotating-axes frame, and it is not displayed in the figure. 



rotating each configuration of a quantum path at imaginary 
time r with respect to the first configuration by i? -1 (r), as 
defined in Eqs. (IT2l . The contribution to the density from 
each quantum path sampled by our RQMC procedure is then 
further rotated, so as to minimize the mutual distance of the 
particle centroids resulting from different paths. In Fig. 3 we 
display the rotating-axes number densiy distribution resulting 
from our simulations for CO@(pH2)i2- 

The results displayed in Figs. 2 and 3 seem to indicate 
that CO@(pH2)i2 clusters behave much as rigid bodies, the 
closest realization of a crystal in a finite system. This behav- 
ior is different from that found in Ref. 9 for undoped (pH2)i3 
clusters whose structure is expected to be very similar to that 
of CO@(pH2)i2, the main difference being the substitution 
of the central CO molecule with another pB.2 molecule. In 
that paper the behavior of (pH2)i3 was described as that of a 
structured superfluid, possibly resembling a supersolid 9 . The 
prediction of a liquid-like behavior for (pH2)i3 — as well as 
its enhanced stability with respect to small clusters of differ- 
ent size — is supported by recent Raman spectroscopy mea- 
surements of cryogenic (pH2)n free jets 13 . Diffusion quantum 
Monte Carlo simulations presented in that paper were used to 
assign individual vibrational lines to clusters of different size, 
and they support the delocalized character of these systems, 
as well as the greatest stability of (pH2)i3- 

In Fig. 4 we report the rotating-axes number-density dis- 
tribution and the I = 6 rotating-axes multipole time correla- 
tion function, resulting from RQMC simulations for (pH2)i3. 
Although the number density distribution displays a nicely 
icosahedral shape — as it is the case for CO@(pH2)i2 — the 
I = 6 multipole correlation function displays a behavior 
which — as it will be shown below — is typical of a molten 
cluster. We interpret this seemingly contradictory behav- 
ior as a consequence of higher frequency of quantum inter- 




FIG. 4: Left panel: pH.2 rotating-axes number density distribu- 
tion for (pH2)i3 (see text). The displayed iso-surface corresponds 
to a value of the density including 60% of the total number of pH 2 
molecules. Right panel: I = 6 rotating-axes multipole time correla- 
tion function (Eq. fT3l for the same system. 



molecular exchanges in (pH2)i3 than in CO@(pH2)i2- This 
is due to both the lack of exchanges between molecules in 
the outer shell and the central molecule in CO@(pH2)i2, and 
also to a larger binding of the molecules in the outer shell to 
the central one, which acts as to hinder intra-shell exchanges. 
During an exchange the shape of the cluster strongly departs 
from icosahedral. We expect therefore that — when the ratio 
between the typical duration of an exchange process and the 
waiting time between two exchanges increases — the I = 6 
multipole autocorrelations would decrease rapidly. It is in- 
teresting to notice that — contrary to CO@(pH2)i2, where the 
decay of the I = 6 autocorrelations is much slower than for 
smaller angular momenta — in the case of (pH2)i3 autocorre- 
lations of different multipoles decay in a qualitatively similar 
way. 




FIG. 5: Left panel: pH.2 rotating-axes number density distribution 
for CO@(pH2)i3 (see text). The displayed iso-surface corresponds 
to a value of the density including 60% of the total number of pFh 
molecules. Right panel: I = 6 rotating-axes multipole time correla- 
tion function (Eq. fT3l for the same system. 

In Fig. 5 we report the rotating-axes number density distri- 
bution and the I = 6 rotating-axes multipole time correlation 
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function, resulting from RQMC simulations for CO@(H2)i3. 
The shape of the number-density iso- surface is strongly irreg- 
ular, as a consequence of the quantum exchanges between one 
of the pH2 molecules in the first solvation shell and the addi- 
tional (thirteenth) molecule at the center of the cluster. The 
I = 6 multipole time correlations correspondigly display a 
rather fast decay. The decay of the these correlations is in fact 
slower than for (£[2)13 for which, however, the shape of the 
number-density isosurface is much more regular (only some- 
what broadened with respect to CO@(H 2 )i2 which we would 
qualify as a solid cluster). 

Conclusions 

We understand that our results — while confirming much of 
what is already known about the quantum behavior of pH 2 
clustera 9 ! 13 ! 18 and their greater propensity to crystallize when 



they are seeded by a foreign molecule 11 — probably raise more 
questions than they answer. A full understanding of the melt- 
ing vs. freezing behavior in these exotic, yet very interest- 
ing, quantum systems will require much more work than it 
has been possible to report in this paper. 
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